Derive growth parameters and natural mortality rates for krill while accounting for spatial heterogeneity in the Western Antarctic Peninsula.
Working Paper to be submitted in a CCAMLR EMM-WG 2024
Mardones, M; Cárdenas, C., Krüger, L., Santa Cruz, F.
24 May, 2024
rm(list = ls())
set.seed(9)
knitr::opts_chunk$set(echo = TRUE,
message = FALSE,
warning = FALSE,
fig.align = 'center',
fig.width = 8,
fig.height = 6,
dev = 'jpeg',
dpi = 300)
#XQuartz is a mess, put this in your onload to default to cairo instead
options(bitmapType = "cairo")
# (https://github.com/tidyverse/ggplot2/issues/2655)
# Lo mapas se hacen mas rapidolibrary(here)
library(kableExtra)
library(ggthemes)
library(ggrepel)
library(ggridges)
#analisis
library(ggpubr)
library(easystats) # multiples unciones analiticas
library(sf)
library(tidyverse, quietly = TRUE)
library(modelsummary)
library(terra) # replace raster
library(TropFishR)
library(mixR)
library(CCAMLRGIS)Methodology
The object ohbio2 come from data exploration analysis in
data request CCAMLR data. This objetc have bio information from
krill.
Join data set with master as c1 set. This join is
trought obs_haul_id variable to get geoposition
variables
ohbio2 <- left_join(c1, ohbio, by="obs_haul_id")
names(ohbio2)## [1] "c1_id.x" "obs_haul_id"
## [3] "obs_logbook_id.x" "obs_haul_number"
## [5] "haul_number.x" "vessel_name"
## [7] "vessel_nationality_code" "fishing_purpose_code"
## [9] "season_ccamlr" "target_species"
## [11] "asd_code" "trawl_technique"
## [13] "catchperiod_code" "date_catchperiod_start"
## [15] "datetime_set_start" "datetime_set_end"
## [17] "datetime_haul_start" "datetime_haul_end"
## [19] "datetime_timezone" "depth_gear_set_end_m"
## [21] "depth_gear_haul_start_m" "depth_bottom_set_end_m"
## [23] "depth_bottom_haul_start_m" "latitude_set_end"
## [25] "longitude_set_end" "latitude_haul_start"
## [27] "longitude_haul_start" "gear_type_code"
## [29] "gear_type" "mesh_code"
## [31] "trawl_net_number" "notes"
## [33] "trawl_duration_depth_h" "trawl_duration_total_h"
## [35] "krill_greenweight_kg" "c1_id.y"
## [37] "obs_logbook_id.y" "haul_number.y"
## [39] "taxon_code" "taxon_scientific_name"
## [41] "taxon_family" "maturity_stage"
## [43] "sex_code" "length_total_cm"
## [45] "greenweight_kg"ohbio3 <- ohbio2 %>%
mutate(Year = year(date_catchperiod_start),
Month = month(date_catchperiod_start),
Day = day(date_catchperiod_start)) %>%
#toupper() para convertir los valores a mayúsculas
mutate(sex_code = toupper(sex_code)) %>%
dplyr::select(7, 9, 11, 12, 14, 24, 25, 29, 42, 44, 46, 47, 43, 45) |>
filter(asd_code=="481")First thing is get different rater layer to join krill data length according different porpoises.
# Cargo linea de costa
coast <- load_Coastline()
coast1<- st_as_sf(coast)
coast2 = st_transform(coast1, "+proj=latlong +ellps=WGS84")
# Uso las agrupaciones de Strata
strata <- st_read("~/DOCAS/Mapas/Antarctic_SHPfiles/Strata.shp",
quiet=T)
strata=st_transform(strata, "+proj=latlong +ellps=WGS84")Test Strata SSMU, just to know another way to join data, but this kind of spatial structuration is deprecated to mamagemente use (Figure@ref(fig:maptest2).
Grouping bio data into stratas
ohbio6 <- st_as_sf(ohbio3 %>%
drop_na(latitude_set_end),
coords = c("longitude_set_end", "latitude_set_end"),
crs = "+proj=latlong +ellps=WGS84")Study Area
# y testeo el mapa
ssmap <- ggplot()+
geom_sf(data = strata |>
filter(ID != "Outer") |>
mutate(ID = str_replace_all(ID, "Extra", "GS")),color="red")+
geom_sf(data = ohbio6 |>
drop_na(length_total_cm) |>
filter(Year>2017),
aes(color = length_total_cm)) +
geom_text_repel(data = strata |>
filter(ID != "Outer") |>
mutate(ID = str_replace_all(ID, "Extra", "GS")),
aes(x = Labx, y = Laby,
label = ID),
min.segment.length = 0,
box.padding = 2,
max.overlaps = 10)+
geom_sf(data = coast2, colour="black", fill=NA)+
scale_color_viridis_c(option = "G",
name="Length (cm)")+
ylim(230000, 2220000)+
xlim(-3095349 , -1858911)+
# coord_sf(crs = 32610)+ #sistema de prpyecccion para campos completos
coord_sf(crs = 6932)+
theme_bw()
ssmap
Strata Maps in 48.1
Length composition by Strata CCAMLR to visualization first. First step is group data into to poligons strata.
strata <- st_make_valid(strata)
sf4 <- st_join(strata, ohbio6)
# Save an object to a file
saveRDS(sf4, file = "sf4.rds")Load RData
# Restore the object
sf4 <- readRDS("sf4.rds")Statistical difference
m <- aov(length_total_cm ~ID, data=sf4)
plot(TukeyHSD(m))
Data
histogram length data to viz in anthor way.
jzstrata <- ggplot(sf4 %>%
mutate(ID = if_else(ID == "Extra", "GERLACHE", ID)) %>%
filter(Year>2010,
ID !="Outer"),
aes(x=length_total_cm,
y = as.factor(Month),
fill=ID))+
geom_density_ridges(stat = "binline", bins = 30,
scale = 1.9,
draw_baseline = FALSE,
alpha=0.9)+
facet_grid(Year~ID) +
geom_vline(xintercept = 3.6, color = "red")+
scale_x_continuous(breaks = seq(from = 1, to = 10,
by = 2))+
scale_y_discrete(breaks = seq(from = 1,
to = 12, by = 4))+
scale_fill_viridis_d(name="Strata",
option="F")+
theme_few()+
xlab("Length (cm.)")+
ylab("")
jzstrata
Now, must filter the DF
sf5 <- sf4 |>
dplyr::select(c("ID",
"date_catchperiod_start",
"length_total_cm",
"Year",
"Month",
"sex_code",
"greenweight_kg")) |>
mutate(ID = if_else(ID == "Extra", "GERLACHE", ID)) |>
filter(ID !="Outer") |>
data.frame()Calculating the proportion of records per year and per stratum in a
table. This table allows us to establish criteria for defining the
strata on which we will analyze krill parameter. We leave out just
JOIN.
# Calcular la proporción
proporcion <- (round(table(sf5$ID) / sum(table(sf5$ID))*100,2))
# Crear una tabla con las proporciones
tabla_proporcion <- as.data.frame(proporcion)
tabla_proporcion$Categoria <- rownames(tabla_proporcion)
rownames(tabla_proporcion) <- NULL
# Renombrar las columnas
colnames(tabla_proporcion) <- c("Strata", "%")
# Mostrar la tabla utilizando kbl()
kbl(tabla_proporcion,
caption = "Proportion length register by Strata") |>
kable_classic(full_width = F,
html_font = "Cambria") |>
kable_styling(bootstrap_options = "striped",
latex_options = "striped")| Strata | % | NA |
|---|---|---|
| BS | 62.71 | 1 |
| EI | 6.58 | 2 |
| GERLACHE | 12.13 | 3 |
| JOIN | 0.14 | 4 |
| SSIW | 18.43 | 5 |
Parameters estimation k and L_inf_
#MALE
sf5fil <- sf5 |>
filter(ID %in% c("BS", "EI", "GERLACHE", "SSIW"),
sex_code %in% "M") |>
mutate(date_catchperiod_start = as.Date(date_catchperiod_start)) |>
mutate(yearly_Group = floor_date(date_catchperiod_start, "year")) |> # New column
drop_na(length_total_cm)
# Definir los nombres correspondientes a cada objeto lfq_results
nombres <- c("SSI", "BS", "GERLACHE", "EI")
# Crear una lista para almacenar los resultados de cada ID
lfq_results <- list()
# Iterar sobre cada ID
for (id in unique(sf5fil$ID)) {
# Obtener el nombre correspondiente
nombre <- nombres[which(unique(sf5fil$ID) == id)]
# Filtrar los datos para el ID actual
df_id <- sf5fil %>% filter(ID == id)
# Crear un objeto lfq para el ID actual
lfq <- lfqCreate(data = df_id,
Lname = "length_total_cm",
Dname = "yearly_Group",
bin_size = 0.1)
# Agregar el resultado a la lista
lfq_results[[id]] <- lfq
# Graficar el objeto lfq
plot(lfq, Fname = "catch",
main = nombre)
}
#FEMALE
sf5filhe <- sf5 |>
filter(ID %in% c("BS", "EI", "GERLACHE", "SSIW"),
sex_code %in% "F") |>
mutate(date_catchperiod_start = as.Date(date_catchperiod_start)) |>
mutate(yearly_Group = floor_date(date_catchperiod_start, "year")) |> # New column
drop_na(length_total_cm)
# Definir los nombres correspondientes a cada objeto lfq_results
nombres <- c("SSI", "BS", "GERLACHE", "EI")
# Crear una lista para almacenar los resultados de cada ID
lfq_resultshe <- list()
# Iterar sobre cada ID
for (id in unique(sf5filhe$ID)) {
# Obtener el nombre correspondiente
nombre <- nombres[which(unique(sf5filhe$ID) == id)]
# Filtrar los datos para el ID actual
df_idhe <- sf5filhe %>% filter(ID == id)
# Crear un objeto lfq para el ID actual
lfqhe <- lfqCreate(data = df_idhe,
Lname = "length_total_cm",
Dname = "yearly_Group",
bin_size = 0.1)
# Agregar el resultado a la lista
lfq_resultshe[[id]] <- lfqhe
# Graficar el objeto lfq
plot(lfqhe, Fname = "catch",
main = nombre)
}
Identified difference in length distribution just to male:
Bhattacharya(lfq_results$EI)## Interactive session needed for Bhattacharya.## NULLBhattacharya(lfq_results$SSIW)## Interactive session needed for Bhattacharya.## NULLBhattacharya(lfq_results$BS)## Interactive session needed for Bhattacharya.## NULLBhattacharya(lfq_results$GERLACHE)## Interactive session needed for Bhattacharya.## NULLnow we assign objet to male lfq_result and female
lfq_resulthe and plot with
lfqRestructure()
Male
# Definir los nombres correspondientes a cada objeto lfq_results
nombres <- c("SSI", "BS", "GERLACHE", "EI")
# Iterar sobre cada objeto lfq almacenado en lfq_results
for (i in seq_along(lfq_results)) {
lfq <- lfq_results[[i]]
# Obtener el nombre correspondiente
nombre <- nombres[i]
# Restructurar el objeto lfq
lfqbin <- lfqRestructure(lfq, MA = 3, addl.sqrt = TRUE)
# Graficar el objeto lfq reestructurado
plot(lfqbin, hist.col = c("white", "black"),
image.col = c(rep(rgb(1,0.8,0.8),1000), "white",
rep(rgb(0.8,0.8,1),1000)),
ylim = c(0,max(lfqbin$midLengths+0.5)),
main = nombre)
# Ajustar curvas al objeto lfq reestructurado
tmp <- lfqFitCurves(lfqbin, par = list(Linf=6.5,
K=0.45,
t_anchor=0.5),
draw = TRUE, col=4, lty=2)
}
Female
# Crear una lista para almacenar los resultados de PW_results <- list()
# Definir los nombres correspondientes a cada objeto lfq_results
nombres <- c("SSI", "BS", "GERLACHE", "EI")
# Iterar sobre cada objeto lfq almacenado en lfq_results
for (i in seq_along(lfq_resultshe)) {
lfqhe <- lfq_resultshe[[i]]
# Obtener el nombre correspondiente
nombre <- nombres[i]
# Restructurar el objeto lfq
lfqbinhe <- lfqRestructure(lfqhe, MA = 3, addl.sqrt = TRUE)
# Graficar el objeto lfq reestructurado
plot(lfqbinhe, hist.col = c("white", "black"),
image.col = c(rep(rgb(1,0.8,0.8),1000), "white",
rep(rgb(0.8,0.8,1),1000)),
ylim = c(0,max(lfqbinhe$midLengths+0.5)),
main = nombre)
# Ajustar curvas al objeto lfq reestructurado
tmp <- lfqFitCurves(lfqbinhe, par = list(Linf=6.5,
K=0.45,
t_anchor=0.5),
draw = TRUE, col=4, lty=2)
}
### Modal Progression analysis with ELEFAN to estimated
Linfand K
Method based on Pauly & Morgan (1987) and Mildenberger et al. (2017)
Krill Male
first we use mixR to identified numbers of compoenent
(modal compositions) we have each stratum. With this result, set
ELEFAN method
unique_ids <- unique(sf5$ID)
obj_id <- list()
for (id in unique_ids) {
filtered_data <- sf5 %>%
drop_na(length_total_cm) %>%
filter(sex_code %in% "M",
ID !="JOIN",
ID == id)
obj_id[[id]] <- filtered_data
}
# select number of components
s_normalgs = select(obj_id$GERLACHE$length_total_cm, ncomp = 2:6)
plot(s_normalgs)
s_normalei = select(obj_id$EI$length_total_cm, ncomp = 2:6)
plot(s_normalei)
s_normalbs = select(obj_id$BS$length_total_cm, ncomp = 2:6)
plot(s_normalbs)
s_normalsswi = select(obj_id$SSWI$length_total_cm, ncomp = 2:6)
plot(s_normalsswi)
# fit a Normal mixture model
modgs = mixfit(obj_id$GERLACHE$length_total_cm, ncomp = 3,
pi = c(0.5, 0.5, 0.5),
mu = c( 3, 4, 5),
sd = c(0.2, 0.2, 0.2))
# fit a Normal mixture model
modei = mixfit(obj_id$EI$length_total_cm, ncomp = 4,
pi = c(0.5, 0.5, 0.5, 0.5),
mu = c( 3, 4, 5, 6),
sd = c(0.2, 0.2, 0.2, 0.2))
# fit a Normal mixture model
modbs = mixfit(obj_id$BS$length_total_cm, ncomp = 5,
pi = c(0.5, 0.5, 0.5, 0.5, 0.5),
mu = c( 2, 3, 4, 5, 6),
sd = c(0.2, 0.2, 0.2, 0.2, 0.2))
# fit a Normal mixture model
modsswi = mixfit(obj_id$SSIW$length_total_cm, ncomp = 4,
pi = c(0.5, 0.5, 0.5, 0.5),
mu = c( 3, 4, 5, 6),
sd = c(0.2, 0.2, 0.2, 0.2))GSplot <- plot(modgs,
theme = "bw",
title ="GERLACHE")
EIplot <- plot(modei,
theme = "bw",
title ="ELEPHANT ISLAND")
BSplot <- plot(modbs,
theme = "bw",
title ="BRANSFIELD STRAIT")
SSWIplot <- plot(modsswi,
theme = "bw",
title ="SSWI")
ggarrange(GSplot,
EIplot,
BSplot,
SSWIplot,
common.legend = TRUE,
ncol=4,
legend="right")#### Parameters to `SSWI`
# run ELEFAN with simulated annealing
res_SAsswi <- ELEFAN_SA(lfq_results$SSIW,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 6,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.001717 maxTime = 30
## Emini is: -0.2570441366
## xmini are:
## 6.418895303 0.7678352515 0.5033732662 0.490536049 0.8187059425
## Totally it used 30.001737 secs
## No. of function call is: 1129# run ELEFAN with genetic algorithm
res_GAsswi <- ELEFAN_GA(lfq_results$SSIW,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 5,
K = 0.5,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
#### Parameters to `EI`
# run ELEFAN with simulated annealing
res_SAei <- ELEFAN_SA(lfq_results$EI,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 6,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.003961 maxTime = 30
## Emini is: -0.2367092272
## xmini are:
## 7.314306291 0.9412915653 0.6182883382 0.7314349413 0.4690019973
## Totally it used 30.00398 secs
## No. of function call is: 1338# run ELEFAN with genetic algorithm
res_GAei <- ELEFAN_GA(lfq_results$EI,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 5,
K = 0.5,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
#### Parameters to `BS`
# run ELEFAN with simulated annealing
res_SAbs <- ELEFAN_SA(lfq_results$BS,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 6,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.008559 maxTime = 30
## Emini is: -0.3182638225
## xmini are:
## 5.593892905 0.604283656 0.1065582037 0.643085584 0.9565868378
## Totally it used 30.00858 secs
## No. of function call is: 1210# run ELEFAN with genetic algorithm
res_GAbs <- ELEFAN_GA(lfq_results$BS,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
#### Parameters to `GERLACHE`
# run ELEFAN with simulated annealing
res_SAgs <- ELEFAN_SA(lfq_results$GERLACHE,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 6,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.001215 maxTime = 30
## Emini is: -0.2869010907
## xmini are:
## 5.789427124 0.7127949446 0.4013525844 0.5469309241 0.1102455743
## Totally it used 30.001234 secs
## No. of function call is: 1600# run ELEFAN with genetic algorithm
res_GAgs <- ELEFAN_GA(lfq_results$GERLACHE,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
#### Krill Female
#### Parameters to `SSWI`
# run ELEFAN with simulated annealing
res_SAsswihe <- ELEFAN_SA(lfq_resultshe$SSIW,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 6,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.016656 maxTime = 30
## Emini is: -0.2316358367
## xmini are:
## 6.003211424 0.6594774183 0.1511300988 0.1896827698 0.9636289105
## Totally it used 30.016675 secs
## No. of function call is: 1157# run ELEFAN with genetic algorithm
res_GAsswihe <- ELEFAN_GA(lfq_resultshe$SSIW,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
#### Parameters to `EI`
# run ELEFAN with simulated annealing
res_SAeihe <- ELEFAN_SA(lfq_resultshe$EI,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 6,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.070971 maxTime = 30
## Emini is: -0.2964366709
## xmini are:
## 5.585338324 0.9580215203 0.2021615843 0.7399474084 0.666007027
## Totally it used 30.07099 secs
## No. of function call is: 1361# run ELEFAN with genetic algorithm
res_GAeihe <- ELEFAN_GA(lfq_resultshe$EI,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
#### Parameters to `BS`
# run ELEFAN with simulated annealing
res_SAbshe <- ELEFAN_SA(lfq_resultshe$BS,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 6,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.011616 maxTime = 30
## Emini is: -0.4349139705
## xmini are:
## 5.923948187 0.8310357623 0.5259462683 0.3320332617 0.281676233
## Totally it used 30.011635 secs
## No. of function call is: 1228# run ELEFAN with genetic algorithm
res_GAbshe <- ELEFAN_GA(lfq_resultshe$BS,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
#### Parameters to `GERLACHE`
# run ELEFAN with simulated annealing
res_SAgshe <- ELEFAN_SA(lfq_resultshe$GERLACHE,
SA_time = 60*0.5,
MA = 5,
agemax = 5,
seasonalised = TRUE, addl.sqrt = TRUE,
init_par = list(Linf = 6,
K = 0.5,
t_anchor = 0.5,
C=0.5,
ts = 0.5),
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1))## Simulated annealing is running.
## This will take approximately 0.5 minutes.
## timeSpan = 30.015046 maxTime = 30
## Emini is: -0.3825402709
## xmini are:
## 5.659959408 0.9045192397 0.9905485124 0.03362630308 0.549246957
## Totally it used 30.015066 secs
## No. of function call is: 1610# run ELEFAN with genetic algorithm
res_GAgshe <- ELEFAN_GA(lfq_resultshe$GERLACHE,
MA = 5,
seasonalised = TRUE,
maxiter = 10,
agemax = 5,
addl.sqrt = TRUE,
low_par = list(Linf = 5,
K = 0.01,
t_anchor = 0,
C = 0,
ts = 0),
up_par = list(Linf = 7.5,
K = 1,
t_anchor = 1,
C = 1, ts = 1),
monitor = FALSE)## Genetic algorithm is running. This might take some time.
Resuming parameters to both sex by strata to krill
#male
GAbs <- unlist(res_GAbs$par)
SAbs <- unlist(res_SAbs$par)
GAei <- unlist(res_GAei$par)
SAei <- unlist(res_SAei$par)
GAsswi <- unlist(res_GAsswi$par)
SAsswi <- unlist(res_SAsswi$par)
GAgs <- unlist(res_GAgs$par)
SAgs <- unlist(res_SAgs$par)
#female
GAbshe <- unlist(res_GAbshe$par)
SAbshe <- unlist(res_SAbshe$par)
GAeihe <- unlist(res_GAeihe$par)
SAeihe <- unlist(res_SAeihe$par)
GAsswihe <- unlist(res_GAsswihe$par)
SAsswihe <- unlist(res_SAsswihe$par)
GAgshe <- unlist(res_GAgshe$par)
SAgshe <- unlist(res_SAgshe$par)
# join
t_k_linf_m <- rbind(GAbs[1:2],
SAbs[1:2],
GAei[1:2],
SAei[1:2],
GAsswi[1:2],
SAsswi[1:2],
GAgs[1:2],
SAgs[1:2])
t_k_linf_f <- rbind(GAbshe[1:2],
SAbshe[1:2],
GAeihe[1:2],
SAeihe[1:2],
GAsswihe[1:2],
SAsswihe[1:2],
GAgshe[1:2],
SAgshe[1:2])
row_names <- c("GA BS", "SA BS", "GA EI",
"SA EI", "GA SSWI", "SA SSWI",
"GA GS", "SA GS")
rownames(t_k_linf_m) <- row_names
rownames(t_k_linf_f) <- row_names
col_names_m <- c("L inf Male", "K Male")
col_names_f <- c("L inf Female", "K Female")
colnames(t_k_linf_m) <- col_names_m
colnames(t_k_linf_f) <- col_names_f
total_para <- round(cbind(t_k_linf_f,
t_k_linf_m),3)
t_k_linf_total <- as_tibble(total_para, rownames = "group")
means <- t_k_linf_total %>%
summarise(across(where(is.numeric), ~ mean(.x, na.rm = TRUE))) %>%
mutate(group = "Mean")
t_k_linf <- bind_rows(t_k_linf_total, means)write_csv(t_k_linf, "parametros_krill.csv")total Table with L_inf_ and K
t_k_linf %>%
separate(group, into = c("METHOD", "STRATA"), sep = " ") %>%
kbl(booktabs = T,
position="ht!",
caption = "Parametres LH in different strata with 2 algoritms to estimation") %>%
kable_styling(latex_options = c("striped",
"condensed"),
full_width = FALSE)| METHOD | STRATA | L inf Female | K Female | L inf Male | K Male |
|---|---|---|---|---|---|
| GA | BS | 6.603000 | 0.385000 | 6.026000 | 0.6860 |
| SA | BS | 5.924000 | 0.831000 | 5.594000 | 0.6040 |
| GA | EI | 5.945000 | 0.616000 | 6.902000 | 0.7990 |
| SA | EI | 5.585000 | 0.958000 | 7.314000 | 0.9410 |
| GA | SSWI | 6.606000 | 0.523000 | 6.470000 | 0.6460 |
| SA | SSWI | 6.003000 | 0.659000 | 6.419000 | 0.7680 |
| GA | GS | 6.011000 | 0.692000 | 6.441000 | 0.8710 |
| SA | GS | 5.660000 | 0.905000 | 5.789000 | 0.7130 |
| Mean | NA | 6.042125 | 0.696125 | 6.369375 | 0.7535 |
Made VB curves by methot and strata by male
separated group
t_k_linf_sep <- t_k_linf %>%
filter(group != "Mean") %>%
separate(group, into = c("METHOD", "STRATA"), sep = " ")# Definir la función de Von Bertalanffy
von_bertalanffy <- function(t, L_inf, K) {
L_inf * (1 - exp(-K * (t-t0)))
}
# Generar una secuencia de tiempos (por ejemplo, de 0 a 5 años)
t0 <- -0.3
time <- seq(0, 7, by = 0.1)
# Calcular longitud esperada para cada zona y tiempo
df_curves_male <- t_k_linf_sep %>%
dplyr::select(c(1:2, 5:6)) %>%
rowwise() %>%
mutate(length = list(von_bertalanffy(time, `L inf Male`, `K Male`))) %>%
unnest(length)
# Agregar el vector de tiempo al DataFrame df_curves
df_curves_male$time <- rep(time, nrow(t_k_linf_sep))
# Graficar con ggplot
male <- ggplot(df_curves_male,
aes(x = time,
y = length,
group=STRATA,
colour=STRATA)) +
geom_line(linewidth=1.1) +
labs(title = "Male",
x = "years",
y = "Length (cm)") +
theme_few()+
facet_wrap(.~METHOD)+
scale_colour_viridis_d(option="F",
name="Strata")Made VB curves by methot and strata by female
# Calcular longitud esperada para cada zona y tiempo
df_curves_female <- t_k_linf_sep %>%
dplyr::select(c(1:4)) %>%
rowwise() %>%
mutate(length = list(von_bertalanffy(time, `L inf Female`, `K Female`))) %>%
unnest(length)
# Agregar el vector de tiempo al DataFrame df_curves
df_curves_female$time <- rep(time, nrow(t_k_linf_sep))
# Graficar con ggplot
female <- ggplot(df_curves_female,
aes(x = time,
y = length,
group=STRATA,
colour=STRATA)) +
geom_line(linewidth=1.1) +
labs(title = "Female",
x = "years",
y = "Length (cm)") +
theme_few()+
facet_wrap(.~METHOD)+
scale_colour_viridis_d(option="F",
name="Strata")both plot
ggarrange(female, male, common.legend = TRUE,
ncol=1,
legend="bottom")
Method Nonlinear Mixed-Effects Models
Fit a nonlinear mixed-effects model (NLMM) to data, via maximum
likelihood using lme4 (Bates et al.
(2015)) to calcularte parametrs
(work in progress)
#We use firts a descomposition methot with library `mixR` and the
sf52020 <- sf5 %>%
filter(sex_code %in% "F",
ID == "EI",
Year =="2020") %>%
drop_na(length_total_cm)
# fit a Normal mixture model
mod1 = mixfit(sf52020$length_total_cm, ncomp = 4,
pi = c(0.5, 0.5, 0.5, 0.5),
mu = c( 3, 4, 5, 6),
sd = c(0.2, 0.2, 0.2, 0.2))
mod2_weibull = mixfit(sf6$length_total_cm, family = 'weibull', ncomp = 3)
s_weibull = select(sf5a$length_total_cm, ncomp = 2:6, family = 'weibull')
s_normal = select(sf5a$length_total_cm, ncomp = 2:6)
plot(s_weibull)
plot(s_normal)
# plot the fitted model# plot the fitted model# plot the fitted model
plot(mod1)
plot(mod2_weibull)
# fit a Normal mixture model (equal variance)
mod1_ev = mixfit(sf6$length_total_cm, ncomp = 2, ev = TRUE)Natural Mortality in Krill
In this exercises different bioanalogic methods are tested with General Algoritm
MALE
MSSWI <- c(res_GAsswi$par, list(agemax = res_GAsswi$agemax))
# use the function M_empirical to estimate natural mortality
Msswi <- M_empirical(Linf = MSSWI$Linf, K_l = MSSWI$K,
tmax = MSSWI$agemax, temp = 5,
tm50 = 2,
method = c("Pauly_Linf", "Hoenig",
"Hoenig",
"AlversonCarney",
"RikhterEfanov"
))
# Bransfield
BS <- c(res_GAbs$par, list(agemax = res_GAbs$agemax))
# use the function M_empirical to estimate natural mortality
Mbs <- M_empirical(Linf = BS$Linf, K_l = BS$K,
tmax = BS$agemax, temp = 5,
tm50 = 2,
method = c("Pauly_Linf", "Hoenig",
"Hoenig",
"AlversonCarney",
"RikhterEfanov"
))
# GERLACHE
GS <- c(res_GAgs$par, list(agemax = res_GAgs$agemax))
# use the function M_empirical to estimate natural mortality
Mgs <- M_empirical(Linf = GS$Linf, K_l = GS$K,
tmax = GS$agemax, temp = 5,
tm50 = 2,
method = c("Pauly_Linf", "Hoenig",
"Hoenig",
"AlversonCarney",
"RikhterEfanov"
))
# Elephand Island
MEI <- c(res_GAei$par, list(agemax = res_GAei$agemax))
# use the function M_empirical to estimate natural mortality
Mei <- M_empirical(Linf = MEI$Linf, K_l = MEI$K,
tmax = MEI$agemax, temp = 5,
tm50 = 2,
method = c("Pauly_Linf", "Hoenig",
"Hoenig",
"AlversonCarney",
"RikhterEfanov"
))FEMALE
MSSWIhe <- c(res_GAsswihe$par, list(agemax = res_GAsswihe$agemax))
# use the function M_empirical to estimate natural mortality
Msswihe <- M_empirical(Linf = MSSWIhe$Linf, K_l = MSSWIhe$K,
tmax = MSSWIhe$agemax, temp = 5,
tm50 = 2,
method = c("Pauly_Linf",
"Hoenig",
"AlversonCarney",
"RikhterEfanov"
))
# Bransfield
BShe <- c(res_GAbshe$par, list(agemax = res_GAbshe$agemax))
# use the function M_empirical to estimate natural mortality
Mbshe <- M_empirical(Linf = BShe$Linf, K_l = BShe$K,
tmax = BShe$agemax, temp = 5,
tm50 = 2,
method = c("Pauly_Linf",
"Hoenig",
"AlversonCarney",
"RikhterEfanov"
))
# GERLACHE
GShe <- c(res_GAgshe$par, list(agemax = res_GAgshe$agemax))
# use the function M_empirical to estimate natural mortality
Mgshe <- M_empirical(Linf = GShe$Linf, K_l = GShe$K,
tmax = GShe$agemax, temp = 5,
tm50 = 2,
method = c("Pauly_Linf",
"Hoenig",
"AlversonCarney",
"RikhterEfanov"
))
# Elephand Island
MEIhe <- c(res_GAeihe$par, list(agemax = res_GAeihe$agemax))
# use the function M_empirical to estimate natural mortality
Meihe <- M_empirical(Linf = MEIhe$Linf, K_l = MEIhe$K,
tmax = MEIhe$agemax, temp = 5,
tm50 = 2,
method = c("Pauly_Linf",
"Hoenig",
"AlversonCarney",
"RikhterEfanov"
))# junto las bases
Total_M_male <- cbind(Mei, Mbs, Mgs, Msswi)
colnames(Total_M_male) <- c("EI", "BS" , "GS", "SSWI")
mean_m <- colMeans(Total_M_male, na.rm = TRUE)
Total_M_Mean_male <- rbind(Total_M_male, Mean = mean_m)
# Agregar los promedios como una nueva fila al final del array
Total_M_femal <- cbind(Meihe, Mbshe, Mgshe, Msswihe)
colnames(Total_M_femal) <- c("EI", "BS" , "GS", "SSWI")
mean_m <- colMeans(Total_M_femal, na.rm = TRUE)
Total_M_Mean_femal <- rbind(Total_M_femal, Mean = mean_m)Total_M_Mean_male %>%
kbl(booktabs = T,
position="ht!",
caption = "Estimated M by Strata for Male") %>%
kable_styling(latex_options = c("striped",
"condensed"),
full_width = FALSE)| EI | BS | GS | SSWI | |
|---|---|---|---|---|
| Alverson and Carney (1975) | 0.6730 | 0.7670 | 0.618 | 0.8030 |
| Hoenig (1983) - Joint Equation | 0.8690 | 0.8690 | 0.869 | 0.8690 |
| Hoenig (1983) - Fish Equation | 0.8470 | 0.8470 | 0.847 | 0.8470 |
| Pauly (1980) - Length Equation | 1.0460 | 0.9830 | 1.128 | 0.9270 |
| Rikhter and Efanov (1976) | 0.7680 | 0.7680 | 0.768 | 0.7680 |
| NA | NA | NA | NA | NA |
| Mean | 0.8406 | 0.8468 | 0.846 | 0.8428 |
Total_M_Mean_femal %>%
kbl(booktabs = T,
position="ht!",
caption = "Estimated M by Strata for Female") %>%
kable_styling(latex_options = c("striped",
"condensed"),
full_width = FALSE)| EI | BS | GS | SSWI | |
|---|---|---|---|---|
| Alverson and Carney (1975) | 0.831 | 1.0710 | 0.762 | 0.9220 |
| Hoenig (1983) - Joint Equation | 0.869 | 0.8690 | 0.869 | 0.8690 |
| Hoenig (1983) - Fish Equation | 0.847 | 0.8470 | 0.847 | 0.8470 |
| Pauly (1980) - Length Equation | 0.920 | 0.6570 | 0.989 | 0.8030 |
| Rikhter and Efanov (1976) | 0.768 | 0.7680 | 0.768 | 0.7680 |
| Mean | 0.847 | 0.8424 | 0.847 | 0.8418 |
Plot to female
# Convertir la matriz a data.frame
df_female <- as.data.frame(Total_M_Mean_femal)
# Agregar la columna de métodos
df_female$Method <- rownames(Total_M_Mean_femal)
# Convertir el data.frame a formato largo
df_long <- df_female %>%
pivot_longer(cols = EI:SSWI,
names_to = "Stratum",
values_to = "Value")
# Dot Plot
m_female<- ggplot(df_long %>%
drop_na() %>%
filter(Method != "Mean"),
aes(x = Stratum,
y = Value,
fill = Method)) +
geom_point(size = 3,
shape = 21,
color = "black") +
geom_text_repel(aes(label = round(Value, 3)),
size = 3,
box.padding = 0.35,
point.padding = 0.5,
segment.color = 'grey50',
min.segment.length = 0,
nudge_y = 0.05,
direction = "both") +
labs(title = "Female",
x = "",
y = "Natural Mortality",
fill = "Method") +
theme_few() +
scale_fill_viridis_d(option = "H") +
ylim(0, 2)Plot to male
# Convertir la matriz a data.frame
df_male <- as.data.frame(Total_M_Mean_male)
# Agregar la columna de métodos
df_male$Method <- rownames(Total_M_Mean_male)
# Convertir el data.frame a formato largo
df_long_male <- df_male %>%
pivot_longer(cols = EI:SSWI,
names_to = "Stratum",
values_to = "Value")
# Dot Plot
m_male<- ggplot(df_long_male %>%
drop_na() %>%
filter(Method != "Mean"),
aes(x = Stratum,
y = Value,
fill = Method)) +
geom_point(size = 3,
shape = 21,
color = "black") +
geom_text_repel(aes(label = round(Value, 3)),
size = 3,
box.padding = 0.35,
point.padding = 0.5,
segment.color = 'grey50',
min.segment.length = 0,
nudge_y = 0.05,
direction = "both") +
labs(title = "Male",
x = "",
y = "Natural Mortality",
fill = "Method") +
theme_few() +
scale_fill_viridis_d(option = "H") +
ylim(0, 2)both plot
ggarrange(m_female, m_male,
common.legend = TRUE,
ncol=1,
legend="right")
Statistical diferences in female
# Realizar la prueba ANOVA
anova_result <- aov(Value ~ Stratum, data = df_long_male %>% drop_na())
summary(anova_result)## Df Sum Sq Mean Sq F value Pr(>F)
## Stratum 3 0.00015 0.00005 0.004 1
## Residuals 20 0.26111 0.01306# Realizar la prueba post hoc de Tukey
tukey_result <- TukeyHSD(anova_result)
# Convertir los resultados de Tukey a un data frame
tukey_df <- as.data.frame(tukey_result$Stratum)
tukey_df$pair <- rownames(tukey_df)
tukey_df <- tukey_df %>%
mutate(Stratum1 = sapply(strsplit(pair, "-"), `[`, 1),
Stratum2 = sapply(strsplit(pair, "-"), `[`, 2))
# Seleccionar las columnas relevantes
tukey_table <- tukey_df %>%
dplyr::select(Stratum1, Stratum2, diff, `p adj`, `lwr`, `upr`) %>%
rename(Difference = diff,
P_value = `p adj`,
Lower_CI = `lwr`,
Upper_CI = `upr`)
kbl(tukey_table,
caption = "Test to differences between strata") |>
kable_classic(full_width = F,
html_font = "Cambria") |>
kable_styling(bootstrap_options = "striped",
latex_options = "striped")| Stratum1 | Stratum2 | Difference | P_value | Lower_CI | Upper_CI | |
|---|---|---|---|---|---|---|
| EI-BS | EI | BS | -0.0062 | 0.9996920 | -0.1908424 | 0.1784424 |
| GS-BS | GS | BS | -0.0008 | 0.9999993 | -0.1854424 | 0.1838424 |
| SSWI-BS | SSWI | BS | -0.0040 | 0.9999171 | -0.1886424 | 0.1806424 |
| GS-EI | GS | EI | 0.0054 | 0.9997963 | -0.1792424 | 0.1900424 |
| SSWI-EI | SSWI | EI | 0.0022 | 0.9999862 | -0.1824424 | 0.1868424 |
| SSWI-GS | SSWI | GS | -0.0032 | 0.9999575 | -0.1878424 | 0.1814424 |
Statistical diferences in female
# Realizar la prueba ANOVA
anova_result_ma <- aov(Value ~ Stratum, data = df_long %>% drop_na())
summary(anova_result_ma)## Df Sum Sq Mean Sq F value Pr(>F)
## Stratum 3 0.00015 0.000048 0.006 0.999
## Residuals 20 0.15347 0.007673# Realizar la prueba post hoc de Tukey
tukey_result_ma <- TukeyHSD(anova_result_ma)
# Convertir los resultados de Tukey a un data frame
tukey_df_ma <- as.data.frame(tukey_result_ma$Stratum)
tukey_df_ma$pair <- rownames(tukey_df_ma)
tukey_df_ma <- tukey_df_ma %>%
mutate(Stratum1 = sapply(strsplit(pair, "-"), `[`, 1),
Stratum2 = sapply(strsplit(pair, "-"), `[`, 2))
# Seleccionar las columnas relevantes
tukey_table_ma <- tukey_df_ma %>%
dplyr::select(Stratum1, Stratum2, diff, `p adj`, `lwr`, `upr`) %>%
rename(Difference = diff,
P_value = `p adj`,
Lower_CI = `lwr`,
Upper_CI = `upr`)
kbl(tukey_table,
caption = "Test to differences between strata") |>
kable_classic(full_width = F,
html_font = "Cambria") |>
kable_styling(bootstrap_options = "striped",
latex_options = "striped")| Stratum1 | Stratum2 | Difference | P_value | Lower_CI | Upper_CI | |
|---|---|---|---|---|---|---|
| EI-BS | EI | BS | -0.0062 | 0.9996920 | -0.1908424 | 0.1784424 |
| GS-BS | GS | BS | -0.0008 | 0.9999993 | -0.1854424 | 0.1838424 |
| SSWI-BS | SSWI | BS | -0.0040 | 0.9999171 | -0.1886424 | 0.1806424 |
| GS-EI | GS | EI | 0.0054 | 0.9997963 | -0.1792424 | 0.1900424 |
| SSWI-EI | SSWI | EI | 0.0022 | 0.9999862 | -0.1824424 | 0.1868424 |
| SSWI-GS | SSWI | GS | -0.0032 | 0.9999575 | -0.1878424 | 0.1814424 |
References
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting
linear mixed-effects models using lme4.
Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01
Mildenberger, T. K., Taylor, M. H., & Wolff, M. (2017). TropFishR: an R package for fisheries analysis with
length-frequency data. Methods in Ecology and Evolution,
8(11), 1520–1527. https://doi.org/10.1111/2041-210X.12791
Pauly, D., & Morgan, G. R. (1987). Length-based methods in fisheries research
(Vol. 13, p. 468).